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ABSTRACT 

Observations of anisotropics in the brightness temperature of the 21 cm line of neutral hydrogen from the period before reionization 
would shed light on the dawn of the first stars and galaxies. In this paper, we use large-scale semi-numerical simulations to analyse 
the imprint on the 21 cm signal of spatial fluctuations in the Lyman-or flux arising from the clustering of the first galaxies. We show 
that an experiment such as the Square Kilometer Array (SKA) can probe this signal at the onset of reionization, giving us important 
information about the UV emission spectra of the first stars and characterizing their host galaxies. SKA-pathfinders with ~ 10% of the 
full collecting area should be capable of making a statistical detection of the 21 cm power spectrum at redshifts z < 20 (corresponding 
to frequencies v > 67 MHz). We then show that the SKA should be able to measure the three dimensional power spectrum as a 
function of the angle with the line of sight and discuss the use of the redshift space distortions as a way to separate out the diff'erent 
components of the 21 cm power spectrum. We demonstrate that, at least on large scales where the hya fluctuations are linear, they 
can be used as a model independent way to extract the power spectra due to these Lya fluctuations. 

Key words. Cosmology: miscellaneous - large-scale structure of Universe - Galaxies: high-redshift - intergalactic medium 



1. Introduction 

The period of formation of the very first stars is one of the 
least understood epochs in the history of the Universe. The 
21 cm spin-flip transition of neutral hydrogen at high redshifts 
has the potential to open a new observational window to study 
this early period, where the most distant, first galaxies reside, 
and even beyond that. Moreover, 21 cm observations may pro- 
vide the only means to detect signatures of the first galaxies at 
z > 15 in the foreseeable future. While the James Webb Space 
Teles cope (JWStQ) will be able to detect galaxies at z < 10 (e.g., 
IStiavelli et al.ll20 Q4), it will not image all of the sources respon- 
sible for the reionization of the Universe. If current estimates in 
the literature are correct, even deep fields with JWST will only 
detect ~ 30% of the sources responsible for reionization or main- 
taining reionization at z ~ 12 ( Salvaterra et al. 2010). At z ~ 20, 
direct emission from stars will be below the threshold of JWST 
observations and new techniques will be required to probe the 
onset of reionization. One possibility is Lyman-CK emitter sur- 
veys over the next decade, but even these surveys will only be 
able to detect the brightest galaxies for selected redshift ranges 
atz^ lQ(e.g.,[Ou chi et al. 2008; Nilsson et al. 2007; Cubv et al] 
2QQ7l:IStark et aLli2007:.Willis et al.,2008:McMahon et al.,.200 



Hibon et aP l2009h . but not the fainter, first galaxies at red- 
shift z > 20. The next generation of large ground based tele- 
scopes, like the European Extremely Large Telescope (E-ELTE|), 
Thirty Meter Telescope (TMTEI), and Giant Magellan Telescope 
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(GMTQ), will only be capable of detecting at most the very 
brightest galaxies at redshifts z ~ 20. Most of the star forma- 
tion however will be in fainter galaxies making it impossible for 
these telescopes alone to survey the full population. In combi- 
nation with 21 cm observations, which are sensitive to the in- 
tegrated Lya emission of the galaxies, it should be possible to 
form a complete census of the first sources. 

With a wide frequency coverage, observations of the 21 cm 
signal will provide a 3 -dimensional tomographic view of the 
inter-galactic medium (IGM), before and during the cosmologi- 
cal reionization (for a review see Furlanetto et al. 2006). As soon 
as the first galaxies appear, at a redshift of z ~ 20 -30 in the stan- 
dard cold dark matter cosmological model (e.g.j Komatsu et af] 
2010), photons emitted at frequencies between the Lya and 
Lyman limit quickly couple the spin temperature of the hy- 
perfine level population to th e IGM gas temper ature t hroug h 
the Wouthuysen-Field eff'ect (IWouthuvsenlll952t lFieldlll959h . 
Provided the IGM has been cooling adiabatically (X-ray heat- 
ing is still negligible during this early stage of galaxy forma- 
tion) , the 21 cm signal will be strong and observed in absorption 
against the cosmic microwave background (CMB). Because the 
Lya coupling depends on the local Lya radiation, which is cor- 
related with galaxies, the 21 cm signal fluctuations will trace the 
location, mass, emission spectra, luminosity and redshift evolu- 
tion of high redshift galaxies. 

Since the first galaxi es represent rare h igh-cr peaks in the 
cosmic density field (e.g. jTrac & Cen|[2007l) . their spatial distri- 
bution have large fluctuations and we expect the Lya contribu- 
tion to the overall 21 cm signal to be appreciable, considerably 
improving the prospects for detection of the Lya fluctuations 
through 21 cm observations. Although the high redshift-end of 
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the reionization process will be inaccessible to the first genera- 
tion experiments, such as the Low Frequency Array (LOFAF0) 
and the Murchison Widefield Array (MWA0), the second gener- 
ation experiments such as the Square Kilometre Array (SKAQ) 
should have enough collecting area to statistically probe the sig- 
nal. 

On the theoretical front, analytical models have been use- 
ful in predicting the pr obable e volution of the high redshift 
signal ( Barkana & Lo^ l2QQ5bl: iPritchard & Furlanettol I2QQ7I: 
iPritchard &Loeb 20d8V but rely on linear approximations and 
require further comparison with simulations. Numerical simu- 
lations, on the other hand, can provide a self-consistent treat- 
ment of the Lya radiative transfer, taking into account eflTects 
such as the scatte ring in the hya line wings (ISemelin et al.ll2QQ7l: 
iBaek et al.ll2009l) . although they are typically slow to run and are 
limited to small volumes (< 100 Mpc/h). Recent developments 
using semi-numerical algorithms have allowed the rapid genera- 
tion of the high reds hift 21 cm signal during the p re-reionization 
epoch (ISantos etal.ll2008l l20ld iMesinger et aLllMoh in large 
boxes, while maintaining the 3-d structure of the signal on large 
scales as seen in the full numerical simulations. 

In this paper we w ill generate fast simulations as described 
in ISantos et"aD (l2010i) to analyze the dependence of the overall 
21 cm signal on several parameters related to the first galaxies 
in the Universe and consider the ability of an experiment like 
SKA to constrain these parameters. The layout of this paper is 
as follows. We start in Section 2 by describing the 21 cm signal 
and the free parameters of the model that aflTect the Lya fluctu- 
ations. In Section 3 we present the experimental setup used for 
SKA, calculating the expected error on the 3-d power spectrum. 
In Section 4 we study the possibility of constraining the signal 
in a model independent way, using the redshift- space distortions 
and the corresponding measurements on the 3-d power spectrum 
as a function of the angle with respect to the line of sight. We 
end with our conclusions and a discussion of the prospects for 
SKA in Sections. 

Throughout this paper where cosmological parameters are 
required we use the standard set of values = 0.28, Q a = 
0.72, = 0.046, if = lOO/zkms-^Mpc"^ (with/z = 0.70), = 
0.96, and (Ts = 0.82, consistent with the latest measurements 
(iKomatsu et aLlHoiQ) . 

2. The 21 cm signal: Lya fluctuations 

The 21 cm brightness temperature corresponds to the change in 
the intensity of the CMB radiation due to absorption or emis- 
sion when it travels through a patch of neutral hydrogen. It is 
given, at an o bserved frequency v in the direction n, by (see e.g. 
ISantos et al.li2008.) 
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where xhi is the fraction of neutral hydrogen (mass weighted), 
dvr/dr is the comoving gradient of the line of sight component 
of the comoving velocity and we use Sa for the fractional value 
of the quantity a (Sa = ^^^) with 6 for the fluctuation in the 
matter density. 
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Size 


Resolution 


Resolution 






(halos) 


{Lya) 


SI 


1000 Mpc 


0.56 Mpc 


1.667 Mpc 


S2 


143 Mpc 


0.09 Mpc 


0.186 Mpc 


S3 (N-body) 


143 Mpc 




0.186 Mpc 



Table 1. Simulations used in the analysis. Middle column shows 
the resolution used to resolve halos in the semi-numerical code, 
while the last column shows the resolution used in the Lya cal- 
culation. 



The spin temperature (Ts) is coupled to the hydrogen gas 
temperature {Tk) through the spin-flip transition, which can 
be excited by collisions or by the absorption of Lya photons 
(Wouthuy sen-Field eff'ect) and we can write: 



Ts 1 + Xtot \ Tk 



(2) 



where Xtot = Xa + Xc is the sum of the radiative and collisional 
coupling parameters and we are already assuming that the color 
temperature of the Lya radiation field at the Lya frequency is 
equal to Tk- Note that although we talk here about Lya radiation, 
in eff'ect we consider the contribution from all the photons up to 
the Lyman limit frequency, since photons redshifting into Lyman 
series res onances can produce Lya photons as a resu lt of atomic 
cascades (lHiratall2006[: IPritchard & Furlanettoll2006h . When the 
coupling to the gas temperature is negligible (e.g. Xtot - ^),Ts ~ 
Ty and there is no signal. On the other hand, for large Xtot, Ts 
simply follows Tk. 

2.1. Simulations 



We use the code SimFast210 described in ISantos et al.l (l2010l) 
to calculate the spin temperature and the corresponding 21 cm 
signal. The code starts by calculating the linear density field in 
a box followed by the corresponding halo and velocity fields. 
Then, it computes the nonlinear density and halo fields, the col- 
lapsed mass distribution, the corresponding star formation rate, 
the IGM gas temperature, the Lya coupling, the collisional cou- 
pling and finally the 21 cm brightness temperature accounting 
for all the contributions including the correction due to red- 
shift space distortions using the nonlinear velocity field. In or- 
der to probe the full range of k space, we generated two end 
to end simulations using the SimFast21 code, one large simula- 
tion with IGpc in size and 1.6 Mpc resolution (SI) and another 
with higher resolution (0.186 Mpc) but smaller 143 Mpc in size 
(simulation S2). This code has the advantage of being modular, 
allowing to replace some of the steps by simulation boxes gen- 
erated by other techniques. Therefore, for comparison, we also 
applied our Lya calcula tion to the output (star formation rate, 
matter density) from thelTrac et al.l (^2008) simulation, thus pro- 
viding a check of the approximation used to generate the halo 
mass function. We call this simulation of the 21 cm signal, S3, 
with a size of 143 Mpc and a resolution of 0.186 Mpc. Table [T] 
shows a summary of the simulations used. All simulations use 
halos with masses down to 10^ M© corresponding to a minimum 
virial temperature of T ~ 10"^ K. Simulation SI cannot resolve 
halos down to these mass scales and we populate the remain- 
ing halos in each (empty) cell using a Poisson sampling biased 
with the underlying density field, giving a mass function con- 
sistent with N-body simulations once non-linear corrections are 
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applied, as described in lSantos et alJ (|2Q1Q|) . With the new ver- 
sion of the code, we also allow a hybrid approach instead of the 
"Poisson" method, where, after the halo finding algorithm, we 
calculate the unresolved collapsed mass in each empty cell using 
the collapsed fr action /con from the extended Press-Schechter 
formula (see e.g. lZahn et al.ll2QQ7b . This is fine as long as we do 
not want to resolve halos within a cell. Note that this is used at 
the cell level, allowing to also apply the non-linear corrections 
to each cell. Both approaches give similar results in terms of the 
21 cm signal, although the later method can be faster as we de- 
crease the halo minimum mass used in the simulation. 

Simulation S3 was run using a hybrid simulation code 
in modeling cosmic reionization (Trac et al. 2008) , incorporat- 
ing N-body, hydrodynamic, and RT algorithms to solve the 
coupled evo lution of the dark matter , bary ons, and radiation 
(iTrac&PenI [20041 l2006l: iTrac &"Cei]l2007 t). It covers a large 
dynamic range and satisfies the requirements of having suffi- 
ciently high resolution to capture small-scale structure, and si- 
mult aneously a suflSciently l arge volume to reduce sample vari- 
ance (iBark ana & Loeb ^2004). The simulation is run in two steps. 
The first step involves running a high-resolution N-body simula- 
tion with 3072^ dark matter particles on an effective mesh with 
11520^ cells in a comoving box, 143Mpc on a side. We iden- 
tified collapsed dark matter halos on the fly using a friends-of- 
friends algorithm, with a linking length b = 0.2 times the mean 
inter-particle spacing, in order to model radiation sources and 
sinks. With a particle mass resolution of 2.68 x I0^h~^ M©, we 
can reliably locate all dark matter halos with virial temperatures 
above the atomic cooling limit (r ~ 10"^ K) with a minimum of 
~ 40 particles ( Heitmann et al.|[2007 ). and half of this collapsed 
mass budget is resolved with > 400 particles per halo. Our halo 
mass functions are in very good agreement w i th other recent 
work (e.g. IReed et all 120071: iLukicet al.l [2Q07t ICohn & Whitd 
l2008h . The second step produces a hydrodynamic -\- RT simula- 
tion with moderate resolution, but incorporating subgrid physics 
modeled using the high-resolution information from the large 
N-body simulation. Radiation sources are prescribed and star 
formatio n rates calculated using the halo model described in 
iTrac & Cen (20 Q7|). We c onsider only Population II stars from 
starbursts ( S chaerei] l2003a) as contributing to the ionizing pho- 
ton budget. The hydro -hRT simulation utilizes equal numbers 
(N = 1536^) of dark matter particles, gas cells, and adaptive 
rays, where for the latter, we track 5 frequencies above the hy- 
drogen ionizing threshold of 13.6 eV. The photo-ionization and 
photo-heating rates for each cell are calculated from the inci- 
dent radiation flux and used in the non-equilibrium solvers for 
the ionization and energy equations. The initial conditions are 
generated with a common white noise field and a linear transfer 
function calculated with CAMB (Lewis et al. 2000). 

In order to calculate the Lya coupling in a given cell, 
the algorithm assumes that the scattering rate is proportional 
to the total Lym an series ffux arriving in that cell (see again 
ISantos et"aD '2010). This Lyman series flux is obtained through 
a 3 -dimensional integration of the comoving photon emissivity, 
6(x, y, z) (defined as the number of photons emitted at position 
X, redshift z and frequency v per comoving volume, per proper 
time and frequency), which is assumed proportional to the star 
formation rate: 



(l2008h simulation. €b(v) is the spectral distribution function of 
the sources (defined as the number of photons per unit frequency 
emitted at v per baryon in stars). We assumed a power law model 
for ebiv): 



e„(x,v,z) = SFRD(x,z)e„(v), 



(3) 



e,(v) = Av", 



(4) 



where SFRD(x, z) is the star formation rate density from the sim- 
ulation (in terms of the number of baryons in stars per comov- 
ing volume and proper time). As stated above, in the simulation 
S3 we used the star formation rate obtained from the lTrac et al.l 



between yQ,(10.2 eV) and the Lyman limit frequency (13.6 eV). 
This will allow us to calculate how sensitive the 21 cm power 
spectrum is to changes in the parameters of the emission model. 

We take a = -0.9 and A such that the integration between the 
Lya and Lyman limit frequencies gives a total emission of 20000 
photons per baryo n. These numbers are based on the expected 
PopII spectra from lSchaereri (2003b). Note that the parameter A 
will be totally degenerate with the star formation rate efl&ciency 
in our model. Since the initial mass function (IMF) is likely to 
be dominated by massive stars, this single power law model is 
likely to be a good d escription of the emission spectra (see e.g. 
iLeitherer et al.lll999b , although a broken power law can some - 
times provide a slightly better fit (iPritchard & Furlanettoll2006h . 
Moreover, as we shall see later, results are very insensitive to 
Of, essentially depending on the total number of photons emitted 
per baryon. We also note that the above calculation of Lya radi- 
ation does not take into account full radiative transfer effects via 
their scattering off of n eutral hydrogen atoms (scattering in the 
wings of the Lya line, Chuzhoy & Zheng 20 071: ISemelin et al.l 
2007, could increase the amplitude of the fluctuations on small 
scales). As a result, we expect that some details of the signal on 
small scales due to Lya coupling to be partly inaccurate. 

Since the details of the galaxy emissivity depend upon the 
assumed initial mass function (IMF), which is poorly known 
( Barkana & Loeb 2001), these values are subject to considerable 
uncertainty. Constraints would therefore provide useful informa- 
tion about the properties of the first galaxies. These uncertainties 
make the mapping between the mean value of the Lya coupling 
(Xa) and redshift uncertain. We will therefore describe the shape 
of the Lya power spectrum for a given (Xa) and caution that the 
redshift at which that actually occurs may be very diflTerent from 
our fiducial model. 



2.2. Comparison 

The 21 cm signal occurs over a wide range of scales, with Lya 
fluctuations contributing at wave-numbers 0.01 < k/(h/Mpc) < 
10. Since none of our simulations can cover this broad set of 
scales individually, we exploit the modularity of the code to 
piece together two simulations to cover the full range. To check 
that this gives self-consistent results, in Figure [T] we plot the full 
21 cm brightness temperature for the three simulations at a few 
values of the mean Lya coupling (Xa). In each case, we take into 
account fluctuations from the density, velocity and Lya terms as 
well as the collisional coupling, ionization fraction (although at 
high redshifts this is negligible) and gas temperature (note again 
that for simulation S3 we calculate the temperature using our 
semi-numerical code based on the simulation SFRD). We see a 
smooth transition from the SI to the S2 simulation, showing for 
the first time, the full range of scales from k ~ 0.01 h/Mpc to 
^ ~ 20 h/Mpc with an amplitude that is roughly between 100 
mK^ and 1000 mK^ on the range of interest, which should help 
the detection of the signal. The increase in the amplitude as the 
coupling increases is due primarily to the cooling of the IGM gas 
in this redshift range, which leads to a stronger absorption signal. 
Note that, as explained in Santos et al. (2010), when integrating 
the Lya flux, we assume the star formation rate is homogeneous 
above a certain large scale in order to make the calculation faster. 



3 



M.G. Santos et al.: Probing the first galaxies with the SKA 




<Xa>= 0.06 

<Xa>= 0.4 

<Xa>= 0.9 









k [h/Mpc] 

Fig. 1. The 21 cm brightness temperature with all fluctuations 
included. Dotted - simulation SI, solid - simulation S2, dashed - 
simulation S3. Redshifts from top to bottom are: 19.25 ((Xa) = 
0.9), 20.25 ({Xa) = 0.4), 22.5 = 0.06). The black dash- 
dotted line shows the correction for (Xa) = 0.4 when doing the 
full Lya flux inhomogeneous integration. 
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Fig. 2. Power spectrum of for several values of (Xq,). Dotted 
lines correspond to simulation SI, solid lines to simulation S2 
and the top dashed line to simulation S3. Redshifts from top to 
bottom are: 19.25 ((jc^) = 0.9), 20.25 «jc^) = 0.4), 22.5 ((jc^) = 
0.06), 18.75 {{Xq) - 8.4). Again, the top, cyan dashed line shows 
the correction when doing the full calculation. 



In this paper this scale was set to 143 Mpc which is the origin of 
the slight break we see in figure [T] for Xq, > 1 . The black dash- 
dotted line in this figure shows the correction for x^ - 0.4 if 
instead we did the inhomogeneous integration up to 300 Mpc. 
Since this is a small correction, we decided to keep the limit to 
143 Mpc. 

For comparison, the dashed line shows the power spectrum 
from simulation S3, which agrees quite well with results using 
the semi-numerical code validating our halo/star formation pre- 
scription (although in S3 there is a decrease of the signal on 
very small scales, which is related to a smoothing of the ve- 
locity term on these scales). An extensive comparison made by 
iBaek et aP feOlO) using radiative transfer simulations (although 
for smaller box size ~ 100 Mpc/h and lower mass resolution 
> lO^^M©) shows a similar spectrum for models in the same 
parameter range. 

Since it will be the focus of the present paper, we now exam- 
ine the Lyof fluctuations in more detail. In Figure |2] we show the 
power spectrum of the quantity = for the three simula- 
tions at four diflTerent values of (Xa). represents the contribu- 
tion of the Lya coupling to the 21 cm signal and varies between 
and 1 . Note that the Xa field tends to be dominated by the large 
values close to individual sources, which obscures the details of 
the coupling since in these regions Aa saturates. The amplitude 
of the power spectrum principally depends on (Xa) with the fluc- 
tuations peaking at (Xq,) ~ 0.5, suggesting that this is the regime 
to best extract information related to the Wouthuy sen-Field ef- 
fect (z ~ 20 in our model), and we will concentrate on these 
redshifts from now on. Figure [2l shows again the correction to 
the calculation if we do not impose a limit to the Lya inhomo- 
geneous integration (top, cyan dashed line, for Xa = 0.9). This 
correction will be more evident as Xa increases (see the "break" 
for Xa = 8.4) but at the same time the contribution to the over- 
all 21cm signal will become smaller so that it should be safe to 
neglect this correction. 

Comparing to figure 14 in lBaek et al.l (l2009l) we see that the 
fluctuation level is similar, although we do not see their small 
scale increase in power due to point sources for k > 1 Mpc/h. 



The diflTerence on small scales (^ > 1 Mpc/h) between our re- 
sults and those of (iBaek et al.ll2009l) is, at least in part, likely due 
to lack of Lya radiative transfer treatment in ours (although dif- 
ferences in size and mass resolution of the simulations should 
also be considered). As will become clear later, it is the 21 cm 
observations at large scales that oflTer the best means to extract 
Lya physics, so the inaccuracies at small scales do not pose a 
significant problem. 

The full N-body simulation (S3) agrees reasonable well with 
our calculation and again we see a smooth transition from the 
very large scale simulation to the high resolution one, giving us 
a complete view over all scales of the expected signal. Although 
it is clear that there is room for more detailed numerical work 
on the 21 cm signal, these comparisons give us confidence in 
the semi- analytic simulation and we now turn to exploring the 
consequences for observations. 

2.3. Identifying the Lya epoch 

In principle, there are several contributions to the 21cm bright- 
ness temperature and full modeling of the signal would be re- 
quired to disentangle the part due to Lya fluctuations and obtain 
parameter constraints on the first galaxies. However, as we can 
see from figure [3] and figure IH both the fluctuations in the col- 
lisional coupling and in the gas temperature should be subdom- 
inant at the redshifts where the Lya fluctuations are strongest 
(when (Xa) < 1 and k < 1/z/Mpc) making the analysis of the 
signal more straightforward (the ionization is very small at these 
redshifts, so that their contribution to the fluctuation power is 
completely negligible). The fluctuations in the gas temperature 
shown in figured only include the heating due to x-rays, which 
can potentially be harder to separate from the Lya fluctuations. 
On top of this we also have fluctuations in the gas temperature 
due to the adiabatic cooling process which are essentially pro- 
portional to the matter density perturbations. Both contributions 
are intrinsically included in the simulation, although, as we see, 
the eflfect from x-ray heating is negligible. 
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Fig. 3. The average and Xc coupling parameters as a function 
of redshift for the 82 simulation. The Lya coupling dominates 
over collisions at all redshifts considered. 



Fig. 5. ^liiSTi) - (STb))^} (rms) of the signal as a function of 
redshift for simulation 82 (frequency in parenthesis). Red lines 
indicate the region where we can safely ignore fluctuations from 
collisions and X-ray heating. 
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Fig. 4. The fluctuations from Lya , Aa (solid) and the gas tem- 

T 

perature, 1 - 7^ (dashed), contributing to the 21 cm brightness 
temperature. Note that both quantities are normalized and the 
gas temperature fluctuations only include the x-ray heating. 

The question then is if whether we can actually identify 
this epoch from observations rather than theoretical guesses. In 
many models, the global (average) 21 cm signal clearly shows 
the transition from Lya coupling to X-ray heating and finally 
to emission once the g as is heated above the CMB temperature 
(iPritchard & Loebl l201Q). This average signal will not be acces- 
sible to interferometers, but as was shown in Santos et al. ( 2008), 
figure 9, the large scale 21 cm power spectrum can also depict 
the evolution of the signal. 

The rms of the signal ({(STb - (STt))^}) can also provide 
that information and should be easier to measure with inter- 
ferometers. As we can see in figure [3 at z > 25 the signal is 
initially small since coupling to the gas temperature is low. As 
the Lya coupling increases, the spin temperature approaches the 
gas temperature and the rms of the signal increases, not only 
due to fluctuations on the Lya field but also because the factor 
Aq, ( 1 - TylTx) is increasing (in absolute terms) since the gas is 



still cooling adiabatically. Once the gas temperature begins to in- 
crease due to X-ray heating the signal decreases and we see the 
turning point at z ~ 16 (the transition from absorption to emis- 
sion occurs later at z ~ 13. Finally, the gas is heated far above 
the CMB temperature and the signal plateaus until reionization 
finally causes it to die away at z ~ 7. 

Note that as X-rays start heating the IGM, the signal does 
not drop immediately since this heating is inhomogeneous and 
there will be some region still cooling adiabatically as well as 
temperature fluctuations contributing to the rms. Therefore we 
see a small plateau at the peak of the signal starting at z ~ 19. 
The conclusion is that it should be safe to neglect the contribu- 
tions of the X-ray heating to the 21 cm signal between the points 
when the rms (or the average) starts rising at z > 24 and before 
we reach the plateau at the maximum z < 19 (Xq, ~ 0.9), at least 
on large scales where the Lya contribution dominates over the 
temperature as we can see in figure IH 

In this case, the equation simplifies and we can write the full 
21 cm signal as: 

dT,{v) = C(z)(l + ^)(1 + SaJH + SaJH +Pd\ (5) 

where (5 = 2/3Tk/(Tk - Ty) (note the diflTerent definition from 
iBarkana & LoebluOOSbh . corresponding to the adiabatic cooling 

process (Tk oc p^'^) and f/^ ^ 180 ({^)^ K. Also, 
X ,,/0.7\/Q,/z2\r/0.15 \ /l+z\l'/' 

(1 - ^){A,){A,) mK 

with Ay = [uTijj^z^rjd?) (^"^v) ~ !)• Figure [6] shows the contribu- 
tion to the full 21 cm signal from each of the terms considered 
on equation[5]using simulations SI and S2. Again we see that the 
Lya term dominates over large scales ^ < 1 h/Mpc, which makes 
it all important to be able to generate large volume simulations 
(> 100 Mpc/h) and give us confidence that it should be possi- 
ble to extract the Lya signal over the redshift range proposed 
without confusion from other contributions. Moreover, the small 
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Fig. 6. The 21 cm power spectrum and its main contributions at 
(Xa) = 0.4 (z ~ 20.25) using simulation SI (red) and simulation 
82 (blue). Top solid lines - all contributions included; dashed - 
Lya only; dotted - matter density; bottom thin solid lines - fluc- 
tuations on the gas temperature from adiabatic cooling only; dot- 
dashed - velocity fluctuations only. 



contribution at these scales from the other terms in equation [5] 
can be modeled with reasonable accuracy if we know the matter 
power spectrum. 

3. Measurements with the SKA 

In this Section, we outline details related to the experimental 
measurement of the 21 cm signal at the very high redshifts dur- 
ing the epoch where Lya dominates. 



3.1. Noise Power Spectrum 

Before we go into the details of the experimental setup, 
we quickly review the noise calculation that goes into our 
analysis. Following iBowman et al.l (l2006h : iMao et al.l (l2008t) : 
iMcQuinn et al.l (l2Q06h . the expected error AP(k, 6) on the mea- 
surement of the 3-d power spectrum of the signal Ps{k, 6) with 
noise Pa^(^, 6) is given by 



AP{K 0) 



1 



VaWM) 



(7) 



where we are assuming that the power spectrum depends only 
on the moduli (k) of the vector k and on the angle 6 between k 
and the line of sight. Nmik, 6) is the total number of modes in k 
space contributing to the measurement (note that the sum is only 
done over half the sphere). In order to obtain the above expres- 
sion, we "grid" the k space into pixels of size dk\dk\\, where _L 
denotes the component perpendicular to the line of sight and || 
the parallel one and assume that measurements in difl'erent pix- 
els are uncorrelated. The size of these pixels is just determined 
by the real space volume of the observed region: 



Vol = ryApovB, 



(8) 



where Apov is the field of view of the experiment (in radi- 
ans), B is the bandwidth for this particular measurement, r{z) = 
£cH~^dz' is the comoving distance to redshift z and y = 



^21 (1 + z)^/H(z) is a conversion factor between frequency in- 
tervals and comoving distances (H(z) is the Hubble expansion 
rate and A21 ~ 21.1cm the rest frame wavelength of the 21 cm 
line). We can relate dk_i_ to the resolution in the w - v (visibility) 
space, du (the Fourier dual of the angular coordinates on the sky) 
by dk_i_ = 2ndu/r (using our Fourier conventions). Note that 



du = 1/ ^|A 



FoV 



(9) 



and is related to the efifective collecting area of one element of 
the interferometer, = A^du^, where A is the wavelength of 
observation (but care must be taken in this relation due to multi- 
beaming). The resolution along the line of sight is just dk\\ = 
27r/(yB) (visibilities are Fourier transformed along the frequency 
coordinate). 

Finally, the noise power spectrum is given by 



du^tik, 0) 



(10) 



where tik, 6) is the time spent observing a given k pixel which is 
just the time spent observing the corresponding pixel of resolu- 
tion du^ on the w - V plane, u = rk^Hln). For a total observation 
time, t{), this will be related to the baseline density distribution, 
^(|u|), through tik,6) = todu^n(\u\) = todu^n (rksm(0) / (In)) as- 
suming the baseline density is rotationally invariant (which is a 
good approximation due to Earth rotation). If we make the fur- 
ther assumption that the baseline density distribution is constant 
on thcu-v plane up to a maximum baseline D^ax (which is not 
the same as assuming an uniform distribution of antennas), we 
have 



^(|u|) = 



' • -"-^ win V 



(11) 



noting that the integration over half the plane should give ap- 
proximately A^^/2 (Na is the number of elements in the interfer- 
ometer). The noise power spectrum then further simplifies to: 



PN(k,0) = r^y- 



(12) 



where Atot is the total collecting area of the telescope. The sys- 
tem temperature was modeled by 



Tsysiz) = 50 + 60 



1+z 
4.73 



2.55 



K 



(13) 



where the first term is the receiver noise temperature (assumed to 
be 50 K) and the second term is the much larger sky temperature 
which is dominated by the Galactic Synchrotron at the redshifts 
of interest. Note in particular that neither the bandwidth or the 
frequency resolution show up in the noise calculation (they can 
be used instead to increase the number of modes available 
to reduce the total error). Moreover, looking at the last expres- 
sion we see that Pn decrease as l/Aj^^ once we fix to and the 
maximum baseline. Once we achieve Ps ~ Pn the only way to 
decrease the total error is through A^^ which does not depend on 
the total collecting area or the baseline distribution. 



3.2. Experimental Setup 

The current generation of radio-interferometers now in construc- 
tion such as the MWA and LOFAR will not have the collecting 
area necessary to probe the high redshift universe (z > 15). This 
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will require a second generation experiment such as the SKA, 
with an observational window between 70 MHz and lOGHz 
and planned to be completed by 2020. Th e SKA will probabiy 
be made up of three different instruments (ISchilizzi et al.ll2007l: 
lFaulknei][20 10): SKA-low, tailored for the low frequency sig- 
nal between 70MHz to 450MHz, probably made up of sparse 
Aperture Arrays so that the collecting area scales as A^; SKA- 
mid between 400MHz and 1.4 GHz, using dense Aperture 
Arrays; and SKA-high, using dishes for frequencies above 1.2 
GHz (note the small overlap between ranges). In this paper, 
we will assume a low frequency "SKA type" experiment with 
a setup capable of probing the high redshift, pre-reionization 
epoch. 

The analysis in the previous Section showed that in order to 
measure the imprint of the first galaxies on the 21 cm signal, 
we need to observe at least up to z = 20, corresponding to a 
frequency of v ~ 68MHz. Allowing for some flexibility, we con- 
sider an experiment that would be able to measure frequencies 
down to 60 MHz. Taking into account the design reference for 
SKA, we assumed a sensitivity of 4000 m^K at 70 MHz (which 
in fact requires a total collecting area of ~ 14 Km^, well above 
the IKm^ used to baptize the telescope) and scaling as around 
these frequencies. Note that the planned collecting area has been 
subject to some change along the years and the current design 
for what is ca lled SKA Phase 1 suggests a sensitivity of 2000 
m^K instead (iGarrett et al.ll2010l) . In order to allow for designs 
with lower sensitivities, we also consider instruments with 20% 
and 10% of the total collecting area of the design reference SKA. 
For the distribution of antennas, the SKADS0 reference design 
assumes they would be collected in 250 stations of 180 m diam- 
eter each, with 66% in a 5 Km diameter core and the rest along 
5 spiral arms out to a 180 Km radius. 

The signal we are looking for dominates on scales ^ < 1 
h/Mpc, requiring baselines up to ~ 10 Km for proper sampling 
at the relevant frequencies, e.g. a resolution of 0.76 arcmin at 
z = 20 (kmax± -1.8 h/Mpc - higher resolution can be achieved 
along the frequency direction although the error would be much 
larger since we would only have modes along the line of sight). 
Taking a reasonable setup, we assumed that 70% of the total 
collecting area would be concentrated within this core of 10 Km 
in diameter. We do not use the rest of the collecting area in the 
noise calculation, assuming instead this would be used for point 
source removal and calibration. The same applies when the total 
collecting area is only 10% and 20% of the reference design. 

For the field of view (FoV) we used 10x10 deg^, which, 
taking z = 20 as reference, sets a resolution of roughly dk_i_ ~ 
4.5 X 10"^h/Mpc. Note that the FoV can be traded with the num- 
ber of beams (in case we want to observe separate fields) and the 
total instantaneous bandwidth (in case we want to probe a larger 
redshift range). Assuming that, in the core, correlations within 
stations are also possible, we considered a minimum baseline of 
50 m giving kmini. ~ 9 x 10"^h/Mpc which is enough to probe the 
scales of interest. The frequency interval for the analysis should 
be chosen carefully at these frequencies since for a fixed fre- 
quency interval, dz increases with redshift. For i nstance, an in- 
terval of 8 MHz used in previous analysis (e.g. iBowman et al.l 
l2007l) corresponds to integrating the signal over Jz ~ 2.5 which 
will average out important features during the period when Lya 
fluctuations are important. On the other hand, using a small in- 
terval means we will have fewer modes along the line of sight. 
Taking a compromise, we assumed that each measurement was 
done with an interval of 4 MHz and a resolution of 0.01 MHz, 
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Parameters — 


Values 


Min. Freq. — 


oU Mhlz 


Sensitivity — 


4(J(JU m /K 


iviaA. uddCiiiic 


10 Km 


Min. baseline — 


50 m 


FoV — 


100 deg2 


Integration time — 


1000 hours 


Freq. interval — 


4 MHz 


Freq. resolution — 


0.01 MHz 



Table 2. Assumed experimental Setup. We also considered 20% 
and 10% of this sensitivity in our analysis. Note that only 70% 
of the assumed collecting area is used in the core of 10 Km di- 
ameter. 



giving dk\\ = 0.08 h/Mpc and kmax\\ = 16.6 h/Mpc. Note that the 
total instantaneous bandwidth should be much larger than this 
(~ 50 MHz), so we can measure several redshifts at the same 
time. Finally, the total integration time was taken to be 1000 
hours (see table [21 for a summary of the experimental setup). 

We followed the error calculation described above for the 
power spectrum assuming a distribution for the antenna/stations 
of the telescope such that a constant baseline density is ob- 
tained. Although this is difficult to achieve in practice since 
the density of stations normally decreases from the center, it 
turns out that this configuration gives the best signal to noise 
in terms of the brightness temperature maps and it should there- 
fore give a good indication of the telescope capabilities to probe 
the 21 cm signal power spectrum. A final note on foregrounds: 
they can be the most damaging factor at these low frequencies 
(together with the calibration issues raised by the ionosphere) al- 
though it is expected that if the foregrounds are smooth enough 
along the frequency direction, they can be simp l y removed by 
fitting out a smooth function (iJelic et al.l l2008l: [Morales et al.l 
2006; Santos et al. 2005). In this analysis we assume that fore- 
ground cleaning was applied on a region of 32 MHz (so as to 
avoid edge eff'ects on our 4 MHz interval) and was success- 
ful on scales smaller than t his 32 MHz, e.g . we ignore scales 
with k < 0.01 h/Mpc (see Harker et al. zOlOj. Calibration issues 
raised by the ionosphere can also be quite damaging in particu- 
lar due to the large field of view of t his lo w frequ ency interfer- 
ometers dCohen & Rott gerin 2^^20091 : iMateiek & MoralesI 120091: 
iLiu et al.ll2010l) and the situation should become more clear once 
we have results from the first generation of EoR experiments. In 
this paper we concentrate on analysing the required minimum 
setup for an SKA type experiment in order to make a statistical 
detection of the signal, neglecting for now the challenges posed 
by the ionosphere. 

3.3. Power Spectrum constraints 

Figure [7] shows the expected error on the 3 dimensional power 
spectrum, P(k) (integrated over 6), assuming the above setup. 
Even if foreground removal aff'ects larger /:-modes than assumed 
here we still have a huge range of measurements available up to 
k - 1 h/Mpc and in particular, measurements are quite good 
on the interval where the Lya contribution is more important. 
On very large scales it should be possible to measure the power 
spectrum even with only 10% of the collecting area. Sample 
variance dominates on large scales while noise (dashed green 
line) is dominating on small scales (the number of modes on 
small scales is limited due to the size of the maximum base- 
lines). Note that we could redistribute the collecting area so that 



7 



M.G. Santos et al.: Probing the first galaxies with the SKA 




105 



k [h/Mpc] 

Fig. 7. Error on the 3d 21cm power spectrum at z=20.25 (Xa = 
0.4) assuming an SKA type experiment. Solid black - signal for 
simulations S1+S2; dashed green - noise power spectrum (equa- 
tion 12) for the full SKA collecting area, 20% and 10% (increas- 
ing amplitude). The expected error taking into account the avail- 
able number of modes (equation 7) is shown with error bars in 
blue, red and yellow respectively. Measurement of the large scale 
power spectrum should be possible as long as the signal ampli- 
tude is not 100 times smaller than current value (this is controlled 
by Xa/(l + Xa)(l - Ty/Ts) so it should be fine for Xa > 0.01). 
On large scales it should be possible to measure the signal even 
with only 10% of the collecting area. 



the noise power spectrum followed more closely the signal. The 
error on large scales is kept small because there are still a large 
number of measurements A^^ due to the high resolution on the 
u - V space from the large field of view (see eq. O. With the 
above measurements of the 3d power spectra, we should be able 
to constrain some of the characteristics of the first galaxies such 
as their emission spectra. Although, as we have seen, the Lya 
signal dominates on larger scales over other contributions, these 
constraints will be particularly strong if we fix the cosmology, 
which should be a reasonable assumption given the expected 
level of constraints from the Planck satellite. 

Up to now we have been assuming a type II emission spec- 
trum with a total of 20, 000 photons per baryon emitted between 
the Lya and Lyman limit frequency. If instead we assume that 
the stars responsible for the Lya coupling are Pop III, we should 
have around 5, 000 photons emitted per baryon used in stars with 
a spectral index ofa = 0.29. Figure[8]shows the results of chang- 
ing the emission model with the dashed curves corresponding to 
the PopIII case (note that due to the lower number of photons, 
these lines also correspond to a lower (Xa) of 0.1. The dotted 
curve corresponds to changing the spectral index to a = 0.29 
while maintaining the total emission to 5,000 photons, which 
will be difl&cult to distinguish under the current noise expecta- 
tions (the total number of emitted photons per baryon is the most 
relevant parameter for the Lya flux). Another possibility is to as- 
sume that Pop III stars are also being formed in halos with mass 
less than 10^ M©. The green dot-dashed line in figure [8] shows 
the eflTect of considering halos down to 10^ M© from simulation 
SI at the same redshift. Note that this has a much higher cou- 
pling and if we consider higher redshifts so that (Xa) ~ 0.1 the 
line will be closer to the dashed one. 
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a = 0.29 Ny = 5000 
a = -0.9 Ny = 5000 
a = 0.9 Ny = 20000 
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Fig. 8. 21 cm temperature power spectrum for a few Lyman 
alpha emission models and simulations SI and S2. Top black 
solid lines correspond to (Xa) = 0.4 while bottom dashed/dotted 
curves to (Xa) = 0.1 (all at z ~ 20.25). The green dot-dashed line 
uses halos down to 10^ M© and a Pop III type emission spectra. 
The error bars in blue corresponds to the full SKA while the ones 
in red to 10% of the collecting area. 



These calculations show that it should be possible for SKA to 
measure (Xa) from observations of the 21 cm power spectrum at 
z ~ 20. Given a detailed model and small error bars this could be 
measured directly from the shape of the power spectrum. More 
generally, if measured over only a small range of wave-numbers, 
the redshift dependence of (Xa) could be extracted by looking at 
the redshift evolution of the a mplitude and slope in a manner 
analogous to that suggested by iLidz et al.l (l2008l) for reioniza- 
tion. Note that, although during Lya domination diflTerent mod- 
els have similar power spectra for the same (Xa) , its evolution 
with redshift could be used to distinguish these models. 

Measurements of (Xa) in a number of diflTerent redshift bins 
would be instrume ntal in constructing a "21 cm Madau plot" 
(iMadau et al.l ll996). The Madau plot shows the evolution of the 
star formation rate as a function of redshift and producing such 
a plot in the very early Universe would be a major accomplish- 
ment of SKA. Measurements of (Xa) do not give this directly, 
since in our model (Xa) is determined by the product of the star- 
formation rate and the emissivity of the galaxies. Breaking this 
degeneracy observationally will be difl&cult and would lead to 
a plot whose absolute normalisation was uncertain, but whose 
shape tracked the star-formation rate (assuming no evolution of 
the g alaxies). However, for many models of early galaxy forma- 
tion (iLeitherer et al.lll99^lBromm & Lars on 2004) the variation 
in the number of Lya photons produced per baryon is only a fac- 
tor of two or so. Thus with a reasonable theoretical prior a useful 
measure of the star-formation history could be established. 



4. Redshift space distortions 

The signal we are trying to probe is actually anisotropic due to 
the redshift space distortions set by the peculiar velocity of the 
HI gas, which is already being taken into account in equation 
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[51 Using a linear approximation, ~ \^fx) \^{x ) ^^^ 



Av ~ 1 - ;^ ^ so to first order the fluctuation in becomes 



1 ^ 
1 + k^a) H dr 



(14) 



The peculiar velocity field originates from inhomogeneities 
in the density field, so that in Fourier space we can use Kaiser's 
approximation ^^(k) = -yu^^(k) (lKaisei]ll987h to write the 
velocity as a function of its dependence in // = cos(^). This linear 
approximation in the velocity contribution allows the brightness 
temperature power spectra to be separated in only three powers 
of yu: 



with 



PAk) = C\z)Ps{k), 



P^2{k) = 2C\z) 
and 



(l+y8W) + 



1 



1 + {Xa) 



(15) 



(16) 



(17) 



PAk) = C^iz) 



(\+pfPs(k) + 



1 



^ 1 + {Xa) 



(18) 



where the power spectra Pab(k) for generic quantities a and b is 
defined as: 



(27rrsUk-W)Pab(k)^ 



{a(k)b(ky} + {b(k)a(ky} 



(19) 



This decomposition can be particularly useful at the high 
redshifts we are considering, where, as can be seen from figure[6j 
the contribution from the velocity compression can be relatively 
strong. By fitting the observed signal to the above polynomial 
(equation [151) at each k, we can in principle separate the cosmo- 
logical and astrophysical contributions to th e brightness temper- 
ature fluctuations (Ba rkana & Loebll2Q05allb ). One can use Pj^i 
and Pj^o to learn about the first sources of radiation while the 

term could be used to measure the matter power spectrum 
directly, with no interference from any other sources of fluctua- 
tions. 



4.1. Assumptions on linearity 

The decomposition discussed in the previous section relies on as- 
sumptions of linearity in both the velocity field and the Lya fluc- 
tuations. On small scales, both of these assumptions may break 
down and we consider both in turn. As we go to small scales, 
fluctuations on the density field increase and the non-linearities 
in the velocity contribution can become important, so that the 
first order approximation to the full velocity term in equation 
[H breaks down. Figure [3 compares the brightness temperature 
power spectrum from the full calculation with the case where 
we assume the first order approximation for the velocity. We see 
that for k > l.O h/Mpc the non-linearities can become important. 
This result implies that for smaller scales the correct brightness 
temperature power spectra has a much more complicated depen- 
dence in jd than the one given by equation [TSl which will compli- 
cate the model i ndependent ex traction of the astrophysical con- 
tributions (.Shaw & Lewisll2QQ8i) . 




k[h/Mpc] 



Fig. 9. Power spectra of brightness temperature fluctuations us- 
ing simulation SI (dotted lines) and using simulation S2 (solid 
lines) for (xa) = 0.4 (z ~ 20.25). Blue lines (top) use the ex- 
act velocity contribution and red lines (bottom) use a first order 
approximation in the velocity contribution. 




k[h/Mpc] 



10° 



Fig. 10. Power spectra of fluctuations in Xa/(l + Xa) built using 
simulation SI for (Xa) = 0.4 (z ~ 20.25). The figure shows the 
the total fluctuation (solid green line) and its first order approxi- 
mation (dashed red line). 



Fluctuations in the hya background too can leave the linear 
regime. The information from the Lya background and conse- 
quently from the physical processes that originated it, is con- 
tained in Sx^ and it can only be directly obtained from the power 
spectrum angular separation in powers of ji for scales where the 
linear approximation of Aq, is valid. To show the limits where the 
first order approximation can be used, we plotted in figure [TOl the 
power spectra of Aa compared to its linear expression, using the 
values from simulation S 1 . The linear approximation should be 
valid for ^< 0.2 h/Mpc. 

Looking at the histogram of Sx^ in figure \TT\ we see that 
there are few points where Sx„ > 1. Although this number is 
small compared to the overall distribution (most cells in the 
simulation have actually Xq, ~ at these high redshifts), they 
are enough to introduce non-linearities on the power spectrum 
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Fig. 11. Detail of the histogram of fluctuations in Sx„ using sim- 
ulation SI for (Xa) = 0.4 (z ~ 20.25). The fluctuations in Sx^ go 
from -50 tifl 350 



calculation. This basically implies that for small scales the 
fluctuations in Xa will have a large dependence on the random 
distribution of the first galaxies so it will only be possible to 
do a linear approximation on Aq, for scales large enough that a 
few small peaks of large intensity do not dominate the power 
spectra of the fluctuation. Some of the strong coupling regions 
coincide with ionized regions which may diminish the efifect of 
strong fluctuations in the brightness temperature power spectra, 
however according with simulations SI, S2 and S3, at the high 
redshifts we are interested (z > 20) the ionized regions are so 
small that including this quantity has no meaningful eflTect in the 
linear approximation used. 

From this study we can conclude that the brightness temper- 
ature power spectrum angular separation in equation [15] should 
be reasonable for ^ < 1 h/Mpc and the linear approximations in 
equations [TT] and [TS] should be valid for larger scales, k < 0.2 
h/Mpc, al lowing the direct extrac tion of information about the 
Lya field (iBarkana & Loebl2005bl) . This is fortunate since it cor- 
responds to the scales where the Lya fluctuations dominate the 
contributions to the brightness temperature (figure [6]) as well as 
the scales where the SKA should perform better (figure O. For 
0.2 < ^ < 1.0 h/Mpc we can still use the angular decomposition 
but equations [T7] and [18] will include higher order terms. Given 
the scales we are interested in, from now on we will concentrate 
on simulation SI (1000 Mpc) which can be used to calculate the 
power spectrum up to ^ ~ 2.3 h/Mpc. 

4.2. Constraints on P(k,ji) 

As we have seen in the previous section, the redshift space dis- 
tortions introduce an anisotropy in the 21 cm power spectrum 
that, in the linear regime, depends on k and even powers of ji. 
More generally non-linearities in the velocity field will intro- 
duce terms with more complicated dependence ( Shaw & Lewis 
|2008|). Measureme nts of P(^, /i) migh t provide more information 
on the Lya signal (iBarkana & Loebl l2()05a) than the spherical 
average P(^) we considered so far. Figure [12] shows the expected 
errors in the measurements of the 3d power spectrum as a func- 
tion of yu, for an SKA type experiment (table O. The calculation 
follows equation |7] where the number of modes in each bin is 




Fig. 12. Expected measurements for the SKA of the 21 cm power 
spectrum as a function of the angle with the line of sight (z = 
20.25, Xa = 0.4). 



now a function of k and yu. We used bins with a logarithmic spac- 
ing in k and linear in ji (with < yu < 1 since the power spec- 
trum only depends on yu^). The number of modes in each bin is 
accounted for by gridding the k space with a resolution along the 
line of sight set by the interval in frequency used for the analysis 
and a resolution perpendicular to the line of sight set by the field 
of view (see section iT2l) . The minimum size of the bins is given 
by this resolution. Note that the number of independent yu bins 
increase with k, although, as can be seen from figure [12] the er- 
rors will also increase due to the experimental noise. Moreover, 
the power spectrum is quite flat for small k which will make it 
harder to extract the polynomial dependence (this is because the 
Pj^Q term is dominating in equation [TSt. 

Although difl&cult as we have seen, we can try to obtain P^o, 
Pj^2 and Pj^4 by fitting equation [15] to the measured P(k, yu) for 
k < 1.0 h/Mpc. This in turn should allow us to obtain direct con- 
straints on a combination of Ps s and Pss • For small k there 
are fewer bins so the uncertainty in P^n will be high even for 
those scales where the errors in P(k,/j.) are small. The expected 
errors on the term can be obtained using a Fisher matrix ap- 
proach (Fisher 1935). For parameters a,b = [P^q,P^2,P^a}, the 
matrix is calculated through 



1 



dPiki^llj) dP(ki,Jlj) 



^ AP(ki,j^j)^ da 



db 



(20) 



where again P(k,jd) is given by equation [T5] and AP(ki,jj.j) is the 
error for each bin. The co variance of the terms is just the 
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Fig. 13. Contributions for the 21 cm temperature power spec- 
trum using simulations (solid lines), built using simulation SI 
for {xa) = 0.4 (z ~ 20.25). Dashed lines: expected errors. Dotted 
lines: expected errors assuming that is known. Solid lines 
from top to bottom are P^o (green), Pj^i (red) and (blue). 

inverse of the Fisher matrix and figure [13] shows the expected 

1-sigma errors ( ^J{F~^\a)^ We see that it should be possible to 
measure with reasonable accuracy both the P^o and P^i terms 
which can already give us interesting constraints on the Lyof 
field. In fact, if we fix the cosmology (e.g. Ps), then we should 
be able to extract information from P^d^^ and P6,„6,„ separately. 
On the other hand, the will be hard to measure, which will 
make the extraction of cosmological information more difl&cult, 
simply because it is an order of magnitude smaller than the other 
terms on the scales we are focusing here. 

In order to improve this last constraint we would have to in- 
crease the collecting area of the instrument and the frequency 
interval used (4 MHz in this case) so to have more modes along 
the line of sight (note however that this will lead to cosmolog- 
ical evolution of the signal along the frequency bin, as already 
discussed). The SKA should be able to measure modes up to 
^ ~ 10 h/Mpc which could have higher values of P^4, however 
for ^ > 1 h/Mpc the angular decomposition is no longer valid 
and the separate measurements of the P(k,/j,) terms will be of 
little use (except maybe for foreground removal). 

In this case, it might be better to just try to fit the parameters 
of the model directly to the averaged power spectrum P(k) us- 
ing the simulations (as shown in figure [8]). If on the other hand 
we assume we already know P^4 with reasonable accuracy then 
the errors on Pj^i will improve considerably (figure [13]) and any 
degeneracies between P^o and Pj^i will be broken. 

4.3. Extracting the Poisson term 

Finally, by measuring the P(k,/j,) terms, we can in principle 
use an estimator to probe the properties of the first sources 
of radiation directly through the Lya field, without having to 
go through a full, model dependent, paramet er fit to the data. 
This is based on the decomposition proposed in^Bark ana & Loebl 
d^OSa), valid in the linear regime. In that case, Sx^ can be ex- 
pressed as a function of its dependence on eflTects correlated and 
uncorrelated with S as: 

SxS^) = W(k)6(k)^Sp(k). (21) 



Fig. 14. Window function for simulation S 1 considering only 
photons emitted below the Lyman beta limit, for (Xa) = 0.4 
(z ~ 20.25). Solid lines are obtained using equation [22] and dot- 
ted lines are obtained using Pss,^ (k) and Ps(k). 

Note that this decomposition is actually quite general and the 
main assumption here is the linearity on S. dp is the Poisson con- 
tribution that arises from the statistical fluctuations in the num- 
ber density of the rare first galaxies (galaxies being discrete ob- 
jects and therefore not tracking the continuous density field per- 
fectly) and which to first order is uncorrelated with the density 
fluctuations. 

The window function, W(k), contains several eflTects by 
which dx„ is related to the underlyi ng matter density fluctua- 
tions, 6 (see iBarkana & Loebll2005bl for details). The strongest 
eflTect between 6 and dx„ is the biasing of the number density of 
galaxies with respect to the density fluctuations by a factor h{z) 
(the average halo bias) so that an overdense region will have a 
factor (1 b{z)6) more sources (Mo & White 1996). This is the 
only term that is implicitly included in our simulations through 
the way we do the 3d integration of the Lya flux. Note how- 
ever that this term compl etely dominates over the other terms in 
IBarkana & Loebl (l2005bh since b{z) ~ 15 at z ~ 20. In that case, 
we can express W(k) as: 

where D(z) is the linear growth function and we are just consid- 
ering photons emitted below the Lyman beta limit for compari- 

SOn, so that Zmaxd) — 

(3 2/27) (l +z)-l (Barkana & Loeb 2005k 
iPritchard & Furlanettoll2QQ7h . 

Using the parametrization of Sx^ we can write PxS^) = 
W^(k)Ps(k) -h Pp(k) and Ps6_,Sk) = W(k)Ps(k), meaning that we 
can isolate W(k) from the simulations (W(k) = Ps5^^(k)/Ps(k)). 
In figure [14] we show the window function obtained from equa- 
tion [22] with values from simul ation SI, which has the sam e 
behavior as the one computed by IPritchard & F urlanettol (l2007h . 
Figure [14] also shows that the window function obtained from 
equation [22] has a higher amplitude for large scales and small 
scales than the one obtained from simulation SI using Ps and 
Pss^^ . On large scales this is probably a consequence of our box 
still not being quite large enough to encompass the rarest, most 
highly biased objects leading to an underestimate of the power 
on large scales. On small scales, this seems to be a consequence 
of the onset of non-linearities in the density and hya fields. On 
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Fig. 15. Power spectra of fluctuations in Xa using simulation S 1 
for {Xa) = 0.4 {z ~ 20.25). Shown are P^^ (solid line), W^Ps 
(dashed line) and (dotted line). 



scales close enough to clusters of sources that the Lyof coupling 
saturates, we might expect to see a steeper fall ofl' in the 21 
cm power as seen here. Nevertheless, W{k) obtained from the 
simulation follows quite closely the expected theoretical values, 
which supports the decomposition given in equation [21] and val- 
idates the signal extraction we discuss next. 

Also, in figure [15] we show P§^^ , W^Ps and Pp, built with Ps, 
Ps and Pss from simulation S 1 . This shows the behavior ex- 
pected from the analytic model of Bark ana & Loebl (|2005b) and 
gives an amplitude of Pp that is broadly consistent with the val- 
ues that we calculate from their model. The Poisson contribution 
dominates for ^ > 0.2 h/Mpc so that there is just a small inter- 
val where we can expect to probe this term before non-linear 
eflTects correlated to the density field become relevant. Note that 
the Poisson contribution flattens on small scales k > l(h/Mpc) 
slightly more than predicted by the analytic model. This is sim- 
ilar to the drop oflT in power on small scales seen in Wik) and 
probably occurs for similar reasons. 

For large enough scales, we can use the observational data to 
obtain the window function through 



Pp.(k) 



= 2x 



1 -hyS-h 



1 



1 -h (Xa) 



W(k) 



(23) 



and the Poisson contribution uncorrected with S: 



Pun-Sik) = Ppo(k) ■ 



4P^) 



C\z) 



1 



1 -h (Xa) 



Pp(k). (24) 



The use of equations [23] and [24] allows us to characterize 
the Lya background and the number and distribution of galax- 
ies at these early redshifts. In figure [161 we used simulation SI 
to show Pun-s and Pj,2 , which are c onsistent with the results ob- 
tained bv lBarkana & Loebl (l2005bl) . Figure [T7] shows the power 
spectra of the Poisson contribution (equation [24]) and the ex- 
pected errors associated with this measurement. As expected, 
due to the high error in P^a, the measurement of this Poisson 
term will be quite difl&cult even with the assumed experimental 
setup (dashed line). If we assume that the matter power spectrum 
is known with reasonable accuracy from CMB data and galaxy 
surveys, then we can combine the measurement of the Pp4 term 
at several k to constrain C^(z) (eq.[T6]), e.g. the mean brightness 



Fig. 16. Power spectra of brightness temperature fluctuations us- 
ing simulation SI for (Xa) = 0.4 (z ~ 20.25). Shown are P^i 
(green dashed line) 2C(z)Ps (blue solid line), and Pun-6 (red dot- 
ted line). 




k[h/Mpc] 



Fig. 17. Power spectra of C^(z) (1/(1 + {x^W Pp(k) (solid line), 
built using simulation SI for (Xa) = 0.4 (z ~ 20.25). Expected 
error (dashed line) and expected error assuming that Ps is known 
(dotted line), so the error in Pp4 comes only from C(z). Dot- 
dashed line assumes that C(z) is known to high accuracy. 



temperature on the sky. Unfortunately the constraint is still large, 
C^(z) ~ (6.42 ± 5.60) X lO^mK^, which translates into large er- 
rors for the Poisson term (dotted lines). Nevertheless, a detection 
will still be possible, with a signal to noise of 2.6 (dashed line) 
and 3.8 (dotted line). 

If we further assume that we combine this result with the 
measurements from a global experiment in order to measure 
the mean brightness temperature with high accuracy, then we 
can extract all the relevant power spectra reasonably well (dot- 
dashed line). 

5. Conclusions 

In this paper, we have made use of a recently developed fast 
semi-numerical code to explore the possibility for 21 cm obser- 
vations of the time of the first galaxies with SKA. We demon- 
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strated that this code allows the simulation of the full range of 
scales likely to be accessible to observations. Here we have fo- 
cused on the impact of 21 cm fluctuations from large fluctuations 
in the Lya coupling, which arise due to the clustering of the rare 
first galaxies and considered diflTerent emission models. 

Measurement of these fluctuations would provide insight 
into the formation of the very first galaxies and would be highly 
complementary with the next generation of large optical/IR tele- 
scopes. We have shown that SKA-pathfinders with ~ 10% of the 
full collecting area should be capable of making a statistical de- 
tection of the 21 cm power spectrum at redshifts z < 20. With 
the full SKA sensitivity this detection would become a measure- 
ment allowing astrophysical properties of the first galaxies to be 
determined. Note that it is crucial that these experiments are able 
to go below 70 MHz in order to probe the signal. 

Observations with an SKA-like instrument would enable the 
determination of (Xa) as a function of redshift from the ampli- 
tude and shape of the 21 cm fluctuations. Although diflTerent 
models with the same (Xa) give similar power spectra, this evolu- 
tion with redshift should allow to distinguish them. Even though 
there is a strong degeneracy between the star formation rate 
and the UV spectral properties of the galaxies themselves, these 
observations would enable constraints to be placed on the star 
formation rate in the earliest generations of galaxies. Although 
crude, the resulting "Madau plot" would be very useful for un- 
derstanding the formation processes of the first galaxies. 

Moving beyond the angle- averaged power spectrum, we 
have investigated the use of redshift- space distortions to separate 
out diflTeren t components of the power spectrum as suggested 
by Barkana & Loebl (l2005bh . These measurements are difl&cult 
due to instrumental limitations and are further complicated by 
non-linearities on small scales. For the first time, we have inves- 
tigated the non-linearity of Lya fluctuations showing that they 
should be important on scales k > 0.2(h/Mpc). Nonetheless, on 
larger scales the detection of these eflTects is possible with SKA 
level sensitivity and should allow to put direct constraints on the 
Lya power spectrum. 

Redshift space distortions further oflfer the possibility of ex- 
tracting components of the power spectrum that do not corre- 
late with the density field. This Poissonian component contains 
direct information about the number density of sources and so 
gives complementary information to the other sources of fluctu- 
ations. We have shown that a detection is possible and if strong 
assumptions about the underlying cosmology are possible and 
combined with information about the mean 21 cm signal, then 
these fluctuations may be picked out clearly by SKA. 

These results illustrate the potential of 21 cm observations to 
shed new light on the astrophysics during the pre-Reionization 
epoch. In the future, our improved knowledge of cosmological 
parameters will provide a firm foundation to pick out the details 
of galaxy formation in the early Universe. SKA and pathfinders 
capable of observing at frequencies v < 100 MHz will begin to 
access this interesting period and transform our understanding 
of the cosmic dawn. 
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